home *** CD-ROM | disk | FTP | other *** search
/ Super Shareware Collection / Super Shareware Collection.iso / os_2 / clisp.zip / FLOATPRI.LSP < prev    next >
Lisp/Scheme  |  1994-02-05  |  22KB  |  420 lines
  1. ;; Ausgabe von Floating-Point-Zahlen mit PRINT und FORMAT
  2. ;; Michael Stoll 10.2.1990 - 26.3.1990
  3. ;; Bruno Haible 8.9.1990 - 10.9.1990
  4.  
  5. ;; Grundgedanken:
  6. ;; Jede Real-Zahl /= 0 repräsentiert ein (offenes) Intervall. Es wird die-
  7. ;; jenige Dezimalzahl mit möglichst wenig Stellen ausgegeben, die in diesem
  8. ;; Intervall liegt.
  9. ;; Um auch große Exponenten zu behandeln, werden Zweier- in Zehnerpotenzen
  10. ;; erst einmal näherungsweise umgerechnet. Nötigenfalls wird die Rechen-
  11. ;; genauigkeit erhöht. Hierbei wird von den Long-Floats beliebiger
  12. ;; Genauigkeit Gebrauch gemacht.
  13.  
  14. (in-package "SYSTEM")
  15.  
  16. ; Stützt sich auf:
  17. ; (sys::log2 digits) liefert ln(2) mit mindestens digits Mantissenbits.
  18. ; (sys::log10 digits) liefert ln(10) mit mindestens digits Mantissenbits.
  19. ; (sys::decimal-string integer) liefert zu einem Integer >0
  20. ;   einen Simple-String mit seiner Dezimaldarstellung.
  21. ; (sys::substring string start [end]) wie subseq, jedoch für Strings schneller.
  22.  
  23. ;; Hauptfunktion zur Umwandlung von Floats ins Dezimalsystem:
  24. ;; Zu einem Float x werden ein Simple-String as und drei Integers k,e,s
  25. ;; berechnet mit folgenden Eigenschaften:
  26. ;; s = sign(x).
  27. ;; Falls x/=0, betrachte |x| statt x. Also oBdA x>0.
  28. ;;   Seien x1 und x2 die nächstkleinere bzw. die nächstgrößere Zahl zu x
  29. ;;   vom selben Floating-Point-Format. Die Zahl x repräsentiert somit das
  30. ;;   offene Intervall von (x+x1)/2 bis (x+x2)/2.
  31. ;;   a ist ein Integer >0, mit genau k Dezimalstellen (k>=1), und es gilt
  32. ;;   (x+x1)/2 < a*10^(-k+e) < (x+x2)/2 .
  33. ;;   Dabei ist k minimal, also a nicht durch 10 teilbar.
  34. ;; Falls x=0: a=0, k=1, e=0.
  35. ;; as ist die Ziffernfolge von a, der Länge k.
  36. (defun decode-float-decimal (x)
  37.   (declare (type float x))
  38.   (multiple-value-bind (binmant binexpo sign) (integer-decode-float x)
  39.     (if (eql binmant 0) ; x=0 ?
  40.       (values "0" 1 0 0) ; a=0, k=1, e=0, s=0
  41.       ; x/=0, also ist sign das Vorzeichen von x und
  42.       ; |x| = 2^binexpo * float(binmant,x) . Ab jetzt oBdA x>0.
  43.       ; Also x = 2^binexpo * float(binmant,x) .
  44.       (let* ((l (integer-length binmant)) ; Anzahl der Bits von binmant
  45.              (2*binmant (ash binmant 1)) ; 2*binmant
  46.              (oben (1+ 2*binmant)) ; obere Intervallgrenze ist
  47.                                    ; (x+x2)/2 = 2^(binexpo-1) * oben
  48.              (unten (1- 2*binmant)) ; untere Intervallgrenze ist
  49.              (untenshift 0)         ; (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
  50.             )
  51.         (when (eql (integer-length unten) l)
  52.           ; Normalerweise integerlength(unten) = 1+integerlength(binmant).
  53.           ; Hier integerlength(unten) = l = integerlength(binmant),
  54.           ; also war binmant eine Zweierpotenz. In diesem Fall ist die
  55.           ; die Toleranz nach oben 1/2 Einheit, aber die Toleranz nach unten
  56.           ; nur 1/4 Einheit: (x+x1)/2 = 2^(binexpo-2) * (4*binmant-1)
  57.           (setq unten (1- (ash 2*binmant 1)) untenshift 1)
  58.         )
  59.         ; Bestimme d (ganz) und a1,a2 (ganz, >0) so, daß
  60.         ; die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
  61.         ; die ganzen a mit a1 <= a <= a2 sind und 0 <= a2-a1 < 20 gilt.
  62.         ; Wandle dazu 2^e := 2^(binexpo-1) ins Dezimalsystem um.
  63.         (let* ((e (- binexpo 1))
  64.                (e-gross (> (abs e) (ash l 1))) ; Ist |e| recht groß, >2*l ?
  65.                g f     ; Hilfsvariablen für den Fall, daß |e| groß ist
  66.                zehn-d  ; Hilfsvariable 10^|d| für den Fall, daß |e| klein ist
  67.                d a1 a2 ; Ergebnisvariablen
  68.               )
  69.           (if e-gross ; Ist |e| recht groß ?
  70.             ; Da 2^e nur näherungsweise gehen kann, braucht man Schutzbits.
  71.             (prog ((h 16)) ; Anzahl der Schutzbits, muß >= 3 sein
  72.               neue-schutzbits
  73.               ; Ziel: 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10.
  74.               (setq g (+ l h)) ; Anzahl der gültigen Bits von f
  75.               ; Schätze d = floor(e*lg(2))
  76.               ; mit Hilfe der Näherungsbrüche von lg(2):
  77.               ; (0 1/3 3/10 28/93 59/196 146/485 643/2136 4004/13301
  78.               ;  8651/28738 12655/42039 21306/70777 76573/254370 97879/325147
  79.               ;  1838395/6107016 1936274/6432163 13456039/44699994
  80.               ;  15392313/51132157 44240665/146964308 59632978/198096465
  81.               ;  103873643/345060773 475127550/1578339557 579001193/1923400330
  82.               ; )
  83.               ; e>=0 : wähle lg(2) < a/b < lg(2) + 1/e,
  84.               ;        dann ist d <= floor(e*a/b) <= d+1 .
  85.               ; e<0  : wähle lg(2) - 1/abs(e) < a/b < lg(2),
  86.               ;        dann ist d <= floor(e*a/b) <= d+1 .
  87.               ; Es ist bekannt, daß abs(e) <= 2^31 + 2^20 .
  88.               ; Unser d sei := floor(e*a/b)-1.
  89.               (setq d (1- (if (minusp e)
  90.                             (if (>= e -970)
  91.                               (floor (* e 3) 10) ; Näherungsbruch 3/10
  92.                               (floor (* e 21306) 70777) ; Näherungsbruch 21306/70777
  93.                             )
  94.                             (if (<= e 22000)
  95.                               (floor (* e 28) 93) ; Näherungsbruch 28/93
  96.                               (floor (* e 12655) 42039) ; Näherungsbruch 12655/42039
  97.               )       )   ) )
  98.               ; Das wahre d wird durch diese Schätzung entweder getroffen
  99.               ; oder um 1 unterschätzt.
  100.               ; Anders ausgedrückt: 0 < e*log(2)-d*log(10) < 2*log(10).
  101.               ; Nun f/2^g als exp(e*log(2)-d*log(10)) berechnen.
  102.               ; Da f < 100*2^g < 2^(g+7), sind g+7 Bits relative Genauigkeit
  103.               ; des Ergebnisses, also g+7 Bits absolute Genauigkeit von
  104.               ; e*log(2)-d*log(10) nötig. Dazu mit l'=integerlength(e)
  105.               ; für log(2): g+7+l' Bits abs. Gen., g+7+l' Bits rel. Gen.,
  106.               ; für log(10): g+7+l' Bits abs. Gen., g+7+l'+2 Bist rel. Gen.
  107.               (let ((f/2^g (let ((gen (+ g (integer-length e) 9))) ; Genauigkeit
  108.                              (exp (- (* e (sys::log2 gen)) (* d (sys::log10 gen))))
  109.                    ))      )
  110.                 ; Das so berechnete f/2^g ist >1, <100.
  111.                 ; Mit 2^g multiplizieren und auf eine ganze Zahl runden:
  112.                 (setq f (round (scale-float f/2^g g))) ; liefert f
  113.               )
  114.               ; Eventuell f und d korrigieren:
  115.               (when (>= f (ash 10 g)) ; f >= 10*2^g ?
  116.                 (setq f (floor f 10) d (+ d 1))
  117.               )
  118.               ; Nun ist 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10 und
  119.               ; f ein Integer ist, der um höchstens 1 vom wahren Wert abweicht:
  120.               ; 10^d * (f-1)/2^g < 2^e < 10^d * (f+1)/2^g
  121.               ; Wir verkleinern nun das offene Intervall
  122.               ; von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
  123.               ; bis (x+x2)/2 = 2^(binexpo-1) * oben
  124.               ; zu einem abgeschlossenen Intervall
  125.               ; von 10^d * (f+1)/2^(g+untenshift) * unten
  126.               ; bis 10^d * (f-1)/2^g * oben
  127.               ; und suchen darin Zahlen der Form 10^d * a mit ganzem a.
  128.               ; Wegen  oben - unten/2^untenshift >= 3/2
  129.               ; und  oben + unten/2^untenshift <= 4*binmant+1 < 2^(l+2) <= 2^(g-1)
  130.               ; ist die Intervall-Länge
  131.               ; = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
  132.               ; = 10^d * ( f * (oben - unten/2^untenshift)
  133.               ;            - (oben + unten/2^untenshift) ) / 2^g
  134.               ; >= 10^d * (2^g * 3/2 - 2^(g-1)) / 2^g
  135.               ; = 10^d * (3/2 - 2^(-1)) = 10^d
  136.               ; und daher gibt es in dem Intervall mindestens eine Zahl
  137.               ; dieser Form.
  138.               ; Die Zahlen der Form 10^d * a in diesem Intervall sind die
  139.               ; mit a1 <= a <= a2, wobei a2 = floor((f-1)*oben/2^g) und
  140.               ; a1 = ceiling((f+1)*unten/2^(g+untenshift))
  141.               ;    = floor(((f+1)*unten-1)/2^(g+untenshift))+1 .
  142.               ; Wir haben eben gesehen, daß a1 <= a2 sein muß.
  143.               (setq a1 (1+ (ash (1- (* (+ f 1) unten)) (- (+ g untenshift)))))
  144.               (setq a2 (ash (* (- f 1) oben) (- g)))
  145.               ; Wir können auch das offene Intervall
  146.               ; von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
  147.               ; bis (x+x2)/2 = 2^(binexpo-1) * oben
  148.               ; in das (abgeschlossene) Intervall
  149.               ; von 10^d * (f-1)/2^(g+untenshift) * unten
  150.               ; bis 10^d * (f+1)/2^g * oben
  151.               ; einschachteln. Hierin sind die Zahlen der Form 10^d * a
  152.               ; die mit a1' <= a <= a2', wobei a1' <= a1 <= a2 <= a2' ist
  153.               ; und sich a1' und a2' analog zu a1 und a2 berechnen.
  154.               ; Da (f-1)*oben/2^g und (f+1)*oben/2^g sich um 2*oben/2^g
  155.               ; < 2^(l+2-g) < 1 unterscheiden, unterscheiden sich a2 und
  156.               ; a2' um höchstens 1.
  157.               ; Ebenso, wenn 'oben' durch 'unten/2^untenshift' ersetzt
  158.               ; wird: a1' und a1 unterscheiden sich um höchstens 1.
  159.               ; Ist nun a1' < a1 oder a2 < a2' , so ist die Zweierpotenz-
  160.               ; Näherung 10^d * f/2^g für 2^e nicht genau genug gewesen,
  161.               ; und man hat das Ganze mit erhöhtem h zu wiederholen.
  162.               ; Ausnahme (da hilft auch keine höhere Genauigkeit):
  163.               ;   Wenn die obere oder untere Intervallgrenze (x+x2)/2 bzw.
  164.               ;   (x+x1)/2 selbst die Gestalt 10^d * a mit ganzem a hat.
  165.               ;   Dies testet man so:
  166.               ;     (x+x2)/2 = 2^e * oben == 10^d * a  mit ganzem a, wenn
  167.               ;     - für e>=0, (dann 0 <= d <= e): 5^d | oben,
  168.               ;     - für e<0, (dann e <= d < 0): 2^(d-e) | oben, was
  169.               ;                nur für d-e=0 der Fall ist.
  170.               ;     (x+x1)/2 = 2^(e-untenshift) * unten == 10^d * a
  171.               ;     mit ganzem a, wenn
  172.               ;     - für e>0, (dann 0 <= d < e): 5^d | unten,
  173.               ;     - für e<=0, (dann e <= d <= 0): 2^(d-e+untenshift) | unten,
  174.               ;                 was nur für d-e+untenshift=0 der Fall ist.
  175.               ; Da wir es jedoch mit großem |e| zu tun haben, kann dieser
  176.               ; Ausnahmefall hier gar nicht eintreten!
  177.               ; Denn im Falle e>=0: Aus e>=2*l und l>=11 folgt
  178.               ;   e >= (l+2)*ln(10)/ln(5) + ln(10)/ln(2),
  179.               ;   d >= e*ln(2)/ln(10)-1 >= (l+2)*ln(2)/ln(5),
  180.               ;   5^d >= 2^(l+2),
  181.               ;   und wegen 0 < unten < 2^(l+2) und 0 < oben < 2^(l+1)
  182.               ;   sind unten und oben nicht durch 5^d teilbar.
  183.               ; Und im Falle e<=0: Aus -e>=2*l und l>=6 folgt
  184.               ;   -e >= (l+2)*ln(10)/ln(5),
  185.               ;   d-e >= e*ln(2)/ln(10)-1-e = (1-ln(2)/ln(10))*(-e)-1
  186.               ;          = (-e)*ln(5)/ln(10)-1 >= l+1,
  187.               ;   2^(d-e) >= 2^(l+1),
  188.               ;   und wegen 0 < unten < 2^(l+1+untenshift) ist unten nicht
  189.               ;   durch 2^(d-e+untenshift) teilbar, und wegen
  190.               ;   0 < oben < 2^(l+1) ist oben nicht durch 2^(d-e) teilbar.
  191.               (when (or (< (1+ (ash (1- (* (- f 1) unten)) (- (+ g untenshift)))) ; a1'
  192.                            a1
  193.                         )
  194.                         (< a2 (ash (* (- f 1) oben) (- g))) ; a2<a2'
  195.                     )
  196.                 (setq h (ash h 1)) ; h verdoppeln
  197.                 (go neue-schutzbits) ; und alles wiederholen
  198.               )
  199.               ; Jetzt ist a1 der kleinste und a2 der größte Wert, der
  200.               ; für a möglich ist.
  201.               ; Wegen  oben - unten/2^untenshift <= 2
  202.               ; ist die obige Intervall-Länge
  203.               ; = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
  204.               ; < 10^d * ((f-1)*oben - (f-1)*unten/2^untenshift) / 2^g
  205.               ; = 10^d * (f-1)/2^g * (oben - unten/2^untenshift)
  206.               ; < 10^d * 10 * 2,
  207.               ; also gibt es höchstens 20 mögliche Werte für a.
  208.             )
  209.             ; |e| ist recht klein -> man kann 2^e und 10^d exakt ausrechnen
  210.             (if (not (minusp e))
  211.               ; e >= 0. Schätze d = floor(e*lg(2)) wie oben.
  212.               ; Es ist e<=2*l<2^21.
  213.               (progn
  214.                 (setq d (if (<= e 22000)
  215.                           (floor (* e 28) 93) ; Näherungsbruch 28/93
  216.                           (floor (* e 4004) 13301) ; Näherungsbruch 4004/13301
  217.                 )       )
  218.                 ; Das wahre d wird durch diese Schätzung entweder getroffen
  219.                 ; oder um 1 überschätzt, aber das können wir leicht feststellen.
  220.                 (setq zehn-d (expt 10 d)) ; zehn-d = 10^d
  221.                 (when (< (ash 1 e) zehn-d) ; falls 2^e < 10^d,
  222.                   (setq d (- d 1) zehn-d (floor zehn-d 10)) ; Schätzung korrigieren
  223.                 )
  224.                 ; Nun ist 10^d <= 2^e < 10^(d+1) und zehn-d = 10^d.
  225.                 ; a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
  226.                 ; a2 sei das größte ganze a < 2^e * oben / 10^d.
  227.                 ; a1 = 1+floor(unten*2^e/(2^untenshift*10^d)),
  228.                 ; a2 = floor((oben*2^e-1)/10^d).
  229.                 (setq a1 (1+ (floor (ash unten e) (ash zehn-d untenshift))))
  230.                 (setq a2 (floor (1- (ash oben e)) zehn-d))
  231.               )
  232.               ; e < 0. Schätze d = floor(e*lg(2)) wie oben.
  233.               ; Es ist |e|<=2*l<2^21.
  234.               (progn
  235.                 (setq d (if (>= e -970)
  236.                           (floor (* e 3) 10) ; Näherungsbruch 3/10
  237.                           (floor (* e 643) 2136) ; Näherungsbruch 643/2136
  238.                 )       )
  239.                 ; Das wahre d wird durch diese Schätzung entweder getroffen
  240.                 ; oder um 1 überschätzt, aber das können wir leicht feststellen.
  241.                 (setq zehn-d (expt 10 (- d))) ; zehn-d = 10^(-d)
  242.                 (when (<= (integer-length zehn-d) (- e)) ; falls 2^e < 10^d,
  243.                   (setq d (- d 1) zehn-d (* zehn-d 10)) ; Schätzung korrigieren
  244.                 )
  245.                 ; Nun ist 10^d <= 2^e < 10^(d+1) und zehn-d = 10^(-d).
  246.                 ; a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
  247.                 ; a2 sei das größte ganze a < 2^e * oben / 10^d.
  248.                 ; a1 = 1+floor(unten*10^(-d)/2^(-e+untenshift)),
  249.                 ; a2 = floor((oben*10^(-d)-1)/2^(-e))
  250.                 (setq a1 (1+ (ash (* unten zehn-d) (- e untenshift))))
  251.                 (setq a2 (ash (1- (* oben zehn-d)) e))
  252.               )
  253.           ) )
  254.           ; Nun sind die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
  255.           ; die ganzen a mit a1 <= a <= a2. Deren gibt es höchstens 20.
  256.           ; Diese werden in drei Schritten auf einen einzigen reduziert:
  257.           ; 1. Enthält der Bereich eine durch 10 teilbare Zahl a ?
  258.           ;    ja -> setze a1:=ceiling(a1/10), a2:=floor(a2/10), d:=d+1.
  259.           ; Danach enthält der Bereich a1 <= a <= a2 höchstens 10
  260.           ; mögliche Werte für a.
  261.           ; 2. Falls jetzt einer der möglichen Werte durch 10 teilbar ist
  262.           ;    (es kann nur noch einen solchen geben),
  263.           ;    wird er gewählt, die anderen vergessen.
  264.           ; 3. Sonst wird unter allen noch möglichen Werten der zu x
  265.           ;    nächstgelegene gewählt.
  266.           (prog ((d-shift nil) ; Flag, ob im 1. Schritt d incrementiert wurde
  267.                  a             ; das ausgewählte a
  268.                 )
  269.             ; 1.
  270.             (let ((b1 (ceiling a1 10))
  271.                   (b2 (floor a2 10)))
  272.               (if (<= b1 b2) ; noch eine durch 10 teilbare Zahl a ?
  273.                 (setq a1 b1 a2 b2 d (+ d 1) d-shift t)
  274.                 (go keine-10-mehr)
  275.             ) )
  276.             ; 2.
  277.             (when (>= (* 10 (setq a (floor a2 10))) a1)
  278.               ; Noch eine durch 10 teilbare Zahl -> durch 10 teilen.
  279.               (setq d (+ d 1)) ; noch d erhöhen, zehn-d wird nicht mehr gebraucht
  280.               ; Nun a in einen Dezimalstring umwandeln
  281.               ; und dann Nullen am Schluß streichen:
  282.               (let* ((as (sys::decimal-string a)) ; Ziffernfolge zu a>0
  283.                      (las (length as)) ; Länge der Ziffernfolge
  284.                      (k las) ; Länge ohne die gestrichenen Nullen am Schluß
  285.                      (ee (+ k d))) ; a * 10^d = a * 10^(-k+ee)
  286.                 (loop
  287.                   (let ((k-1 (- k 1)))
  288.                     (unless (eql (schar as k-1) #\0) (return)) ; eine 0 am Schluß?
  289.                     ; ja -> a := a / 10 (wird aber nicht mehr gebraucht),
  290.                     ; d := d+1 (wird aber nicht mehr gebraucht),
  291.                     (setq k k-1) ; k := k-1.
  292.                 ) )
  293.                 (when (< k las) (setq as (sys::substring as 0 k))) ; Teilstring nehmen
  294.                 (return (values as k ee sign))
  295.             ) )
  296.             ; 3.
  297.             keine-10-mehr
  298.             (setq a
  299.               (if (eql a1 a2)
  300.                 ; a1=a2 -> keine Frage der Auswahl mehr:
  301.                 a1
  302.                 ; a1<a2 -> zu x nächstgelegenes 10^d * a wählen:
  303.                 (if e-gross
  304.                   ; a = round(f*2*binmant/2^g/(1oder10)) (beliebige Rundung)
  305.                   ;   = ceiling(floor(f*2*binmant/(1oder10)/2^(g-1))/2) wählen:
  306.                   (ash (1+ (ash
  307.                              (let ((z (* f 2*binmant)))
  308.                                (if d-shift (floor z 10) z)
  309.                              )
  310.                              (- 1 g)
  311.                        )   )
  312.                        -1
  313.                   )
  314.                   ; |e| klein -> analog wie oben a2 berechnet wurde
  315.                   (if (not (minusp e))
  316.                     ; e>=0: a = round(2^e*2*binmant/10^d)
  317.                     (round (ash 2*binmant e)
  318.                            (if d-shift (* 10 zehn-d) zehn-d) ; 10^d je nach d-shift
  319.                     )
  320.                     ; e<0, also war d<0, jetzt (wegen Schritt 1) d<=0.
  321.                     ; a = round(2*binmant*10^(-d)/2^(-e))
  322.                     (ash (1+ (ash
  323.                                (* 2*binmant
  324.                                   (if d-shift (floor zehn-d 10) zehn-d) ; 10^(-d) je nach d-shift
  325.                                )
  326.                                (+ e 1)
  327.                          )   )
  328.                          -1
  329.                     )
  330.             ) ) ) )
  331.             (let* ((as (sys::decimal-string a)) ; Ziffernfolge von a
  332.                    (k (length as)))
  333.               (return (values as k (+ k d) sign))
  334. ) ) ) ) ) ) )
  335.  
  336. ; Ausgabefunktion für PRINT/WRITE von Floats:
  337. (defun write-float (stream arg #| &optional (plus-sign-flag nil) |# )
  338.   (unless (floatp arg)
  339.     (error #+DEUTSCH "Argument ist kein Float: ~S"
  340.            #+ENGLISH "argument is not a float: ~S"
  341.            #+FRANCAIS "L'argument n'est pas un FLOAT : ~S"
  342.            arg
  343.   ) )
  344.   (multiple-value-bind (mantstring mantlen expo sign)
  345.       (decode-float-decimal arg)
  346.     ; arg in Dezimaldarstellung: +/- 0.mant * 10^expo, wobei
  347.     ;  mant die Mantisse: als Simple-String mantstring mit Länge mantlen,
  348.     ;  expo der Dezimal-Exponent,
  349.     ;  sign das Vorzeichen (-1 oder 0 oder 1).
  350.     (if (eql sign -1) ; arg < 0 ?
  351.       (write-char #\- stream)
  352.       #| (if plus-sign-flag (write-char #\+ stream)) |#
  353.     )
  354.     (let ((flag (<= -2 expo 7))) ; arg=0 oder 10^-3 <= (abs arg) < 10^7 ?
  355.       ; Was ist auszugeben? Fallunterscheidung:
  356.       ; flag gesetzt -> "fixed-point notation":
  357.       ;   expo <= 0 -> Null, Punkt, -expo Nullen, alle Ziffern
  358.       ;   0 < expo < mantlen ->
  359.       ;     die ersten expo Ziffern, Punkt, die restlichen Ziffern
  360.       ;   expo >= mantlen -> alle Ziffern, expo-mantlen Nullen, Punkt, Null
  361.       ;   Nach Möglichkeit kein Exponent; wenn nötig, Exponent 0.
  362.       ; flag gelöscht -> "scientific notation":
  363.       ;   erste Ziffer, Punkt, die restlichen Ziffern, bei mantlen=1 eine Null
  364.       ;   Exponent.
  365.       (if (and flag (not (plusp expo)))
  366.         ; "fixed-point notation" mit expo <= 0
  367.         ; erst Null und Punkt, dann -expo Nullen, dann alle Ziffern
  368.         (progn
  369.           (write-char #\0 stream)
  370.           (write-char #\. stream)
  371.           (do ((i expo (+ i 1)))
  372.               ((eql i 0))
  373.             (write-char #\0 stream)
  374.           )
  375.           (write-string mantstring stream)
  376.           (setq expo 0) ; auszugebender Exponent ist 0
  377.         )
  378.         ; "fixed-point notation" mit expo > 0 oder "scientific notation"
  379.         (let ((scale (if flag expo 1)))
  380.           ; Der Dezimalpunkt wird um scale Stellen nach rechts geschoben,
  381.           ; d.h. es gibt scale Vorkommastellen. scale > 0.
  382.           (if (< scale mantlen)
  383.             ; erst scale Ziffern, dann Punkt, dann restliche Ziffern:
  384.             (progn
  385.               (write-string mantstring stream :end scale)
  386.               (write-char #\. stream)
  387.               (write-string mantstring stream :start scale)
  388.             )
  389.             ; scale>=mantlen -> es bleibt nichts für die Nachkommastellen.
  390.             ; alle Ziffern, dann scale-mantlen Nullen, dann Punkt und Null
  391.             (progn
  392.               (write-string mantstring stream)
  393.               (do ((i mantlen (+ i 1)))
  394.                   ((eql i scale))
  395.                 (write-char #\0 stream)
  396.               )
  397.               (write-char #\. stream)
  398.               (write-char #\0 stream)
  399.           ) )
  400.           (decf expo scale) ; der auszugebende Exponent ist um scale kleiner.
  401.       ) )
  402.       ; Nun geht's zum Exponenten:
  403.       (let ((e-marker
  404.               (cond ((case *READ-DEFAULT-FLOAT-FORMAT*
  405.                        (SHORT-FLOAT (short-float-p arg))
  406.                        (SINGLE-FLOAT (single-float-p arg))
  407.                        (DOUBLE-FLOAT (double-float-p arg))
  408.                        (LONG-FLOAT (long-float-p arg))
  409.                      ) #\E )
  410.                     ((short-float-p arg) #\s)
  411.                     ((single-float-p arg) #\f)
  412.                     ((double-float-p arg) #\d)
  413.                     ((long-float-p arg) #\L)
  414.            )) )
  415.         (unless (and flag (eql e-marker #\E)) ; evtl. Exponent ganz weglassen
  416.           (write-char e-marker stream)
  417.           (write expo :base 10 :radix nil :stream stream)
  418. ) ) ) ) )
  419.  
  420.